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Abstract 

We have dramatically extended the zero field susceptibility series at both high and 
low temperature of the Ising model on the triangular and honeycomb lattices, and 
used these data and newly available further terms for the square lattice to calculate a 
number of terms in the scaling function expansion around both the ferromagnetic and, 
for the square and honeycomb lattices, the antiferromagnetic critical point. 

Cyril Domb was a pioneer in the application of series expansions to the study of critical 
phenomena [U [2] . He encouraged many colleagues to develop this approach and headed a 
group, the "Kings College group," who applied his ideas to investigate the behaviour of co- 
operative assemblies and percolation processes with considerable success. Bomb's unselfish 
and generous attitude in urging people to follow up and develop the series approach was an 
important factor in the subsequent evolution of research in these areas. It is therefore with 
considerable pleasure that we dedicate this paper to Cyril Domb, on the occasion of his 
90th birthday. In it we show just how powerful the series approach can be, as we present 
an analysis based on hundreds, and in some cases thousands of terms in the expansion 
of the susceptibility of the two-dimensional Ising model. It would be fair to say that no 
method other than the series method provides anything remotely approaching this level of 
information about the susceptibility. 

1 Introduction 

A decade ago a number of the current authors reported on a substantial extension of the 
square lattice Ising susceptibility series to some 300 terms O Hj . We found breakdown of 
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the simple scaling picture that assumes the absence of irrelevant scaling fields. The first 
breakdown, which was identified with the breakdown of rotational symmetry of the square 
lattice, occurred at O(r^), with r to leading order proportional to the temperature deviation 
from critical, T—Tc- A second breakdown was identified at O(t^), ascribed to an additional 
irrelevant variable. At the time it was foreshadowed that the corresponding calculation for 
the triangular and honeycomb lattices would be necessary in order to distinguish between 
lattice effects and more fundamental breakdowns intrinsic to the model. 

In this study we report on the derivation and analysis of triangular and honeycomb 
lattice series to more than 300 terms, followed by a calculation of the corresponding scal- 
ing functions. Our numerical work is of sufficient accuracy that we can unambiguously 
identify the same irrational constant, that appeared at O(r^) in the square lattice scaling 
function and was ascribed to a second irrelevant variable, as a contribution to O(r^) in 
both triangular and honeycomb lattices. Furthermore, we find another irrational constant 
common to all lattices at 0(r^'^) which can be ascribed to yet another (third) irrelevant 
variable. These results clearly indicate aspects of universality in the susceptibility beyond 
those found at leading order. 

A limited selection of our results which are the basis for these remarks on universality 
are given in the immediately following text while the very extensive complete listing can 
be found in the appendices. In subsequent sections we elaborate on the results below and 
give details of how they were obtained. Specifically, in section [2] we put our results in the 
context of scaling theory and speculate on the identification of our correction to scaling 
terms with the operators of the conformal field theory that describes the Ising model. 
Section [3] describes how the series expansions were obtained from the quadratic recurrence 
relations for the Z-invariant Ising model specialized to the triangular /honeycomb system. 
In section |4] we describe some of the series analysis details, in particular those aspects that 
differ from what was done in [3]. 

Our numerical work indicates that the reduced susceptibility on any lattice near the 
ferromagnetic critical point (for T > Tc or T < Tc) is given b}|^ 

-lattice 1 rp lattice /^lattice\ \ — 7 / A rplattice , Tjlattice (-\\ 

where B is the contribution of the "short-distance" terms and includes an analytic back- 
ground. It is of the form 

00 Lv^J 

i? = ^^6fe^)(log|T|rr'' (2) 

q=0 p=0 

with the the same above and below Tc but, of course, different for each lattice. The 
temperature variable r is simply related to the low-temperature elliptic parameter k {= k^) 

^The notation here differs from tliat in [3] and tlie earlier literature in tfiat for a common treatment of 
all lattices it is convenient to absorb a factor (2Kc\^y^'^ into the definition of Co±, cf. equation ^ vs. the 
appendix in [3]. 
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by the same expressions 



^''Vk-^), k = {T+Vl + T^f (3) 



2 V Vk 



for every lattice. The elUptic parameter k depends on the lattice; we have, with K 
J/kBT, 

ksq = l/s^, s = sinh2i^sq, square, 

4^3/2 



exp(— 4i^tr)> triangular. 



khc = —7z TTTT— — ^, z = exp (-2Khc), honeycomb. (4) 

(1 - zy{i + z) 

Duality relates the high-temperature elliptic parameter A;> to the low-temperature one 
by ky = 1/A;< or, what is equivalent, by the replacement r — t- — r. Furthermore, since 
the honeycomb lattice is the dual of the triangular lattice and is also related by a star- 
triangle transformation, we can take ktr = k^c as a common elliptic parameter A:< with the 
u (triangle) and z (honeycomb) then connected by 

~ 1 - z + z2 ' i + ^+^(i_^)(i + 3^- 

The Co± constants in ([T]) for the different lattices are related as follows. First, we define 
Co± as the values for the square lattice, that i^ 

Co+ = Cq\ = 1.00081526044021264711947636304721023693753492559778\ 

92751083189882604491051665192385157187485052515870678 V2, (6) 

Co_ = Co'i = 1.0009603287252621894809349551720973205725059517701173\ 

61531948595158755619871466228353934981038826872108 V2/{l2ir). 

Then 

C*^± = 4Co±/%/27, Co^l = 8Co±/^/27, (7) 
as follows from lattice-lattice scaling |1H [12] or Z-invariance [13] . The scaling functions 



^For the calculation of Co± see the footnote on page 3904 in [5]. Here we have used predictor-correctors 
of order as high as 25. This approach uses the Painleve III equation of [6l[71[8]. Alternatively, one can also 
use the Painleve V formulation El [101 ■ 
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through 0(t ) are 



where 



r2 ( 647 
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(8) 



Ce- = 4.54530659737804996885745146127924976519048127125911619\ 

2274173103880744339809, 
C6+ = 0.118322588863244285519212856456397718968975725227410541191067925, 
Cio- = 0.464207706785944087396503330097938832697360392193891710489569762, 
Cio+ = 0.0123440983021588166317669811773152519959150566201343. (9) 

We have not yet been able to identify these constants but expect them to be of a similar 
status to the constants Co± in ([6]) which are related to solutions of the Painleve III O [71 [8] 
or Painleve V equation |9l [10]. We note that the constants must relate to the expansion 
coefficients in (2.27) of [H], which have to satisfy a Painleve V hierarchy of differential 
equations and should lead to further coefficients Ci2,±, Ci4^-|-, etc. We also note that in Q 
we have split off a factor k^l'^^ leaving only even powers of r in the expansions of F jk}-^. 

The staggered susceptibility at the ferromagnetic point of a bipartite lattice, or what is 
equivalent, the susceptibility for an antiferromagnet, is given by an expression of the same 
form as ([l|. For the square lattice the i^i]'^^ vanishes; there is only a background B^^ as 
found in [1]. On the other hand the Fisher |15| relation 



he 



together with ([T]) and (|8]) implies that if we define 

8Co±/\^ 



Che laf 
Oil 



then 



^he I af 



[-Ce±T^ + Ce±T^ - (C76± - Cio±y + 0(t 



12^ 



(10) 



(11) 



(12) 
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Also, 



B 



hciaf 



2^tr _ ^hc_ (-^3) 

Equations (12) and (13) have provided significant tests confirming the correctness and 



accuracy of our numerical analyses. 

To the constants Cq± in ([9]) one could add the same rational above and below Tc and, 
on absorbing this change in other rationals in ([8|, leave those equations unchanged in form. 
A corresponding replacement Cio± — Cio± + rational x Cq±+ rational with similar con- 
sequences is possible. This non-uniqueness in ([sj ) has been removed by arbitrarily adopting 
the particularly simple form for i^^cjaf (12^7 Note however that any such redefinitions 
can never eliminate the irrationals from ([8) or (12) and we conclude that this is evidence 
for at least two irrelevant scaling fields beyond the one breaking rotational invariance and 
contributing first at O(t^) to the square lattice susceptibility. Furthermore, the presence 
of the same irrationals in the scaling functions in ([s]) is evidence for a universality in terms 
beyond the leading order. We will elaborate on this in section [2] where, among other things, 
we make comparisons with the Aharony and Fisher [16\ [T7] scaling functions. 

It is also possible, based on existing results, to derive the reduced susceptibility of the 
Ising model on the kagome lattice. This is given by [18| eqn. 2.1], in terms of the reduced 
susceptibility of the model on the honeycomb lattice. Further aspects of this connection 
can be found in [19j. With Q = J/k^T for the kagome lattice and z = 2/(e^<3 + 1), this 
equation can be written as 

^ka ^ _ ^2) -he + 1 ((1 + _ (1 _ z^){a,a,)t) . (14) 

We note that the z variable is the same variable as in ([4]), pertaining to the interaction 
strength on the honeycomb lattice that results from reversing the star-triangle and decora- 
tion transformations on the kagome lattice. We can also associate with the kagome lattice 
the elliptic parameter and temperature variable associated with the honeycomb lattice, as 
given in ^ and The average (crjcrj)^^ in (14) is the nearest-neighbour correlation 



function of the honeycomb lattice, which is a simple multiple of the internal energy. It is 



given explicitly by eq. (27) below. 



As the second term in (14) is a "short-distance" term, it does not contribute to the 



scaling function F^, which is thus entirely derived from the first term in (14). By absorbing 
an extra normalising factor associated with 1 — z'^ into the constant term, we derive 



r' 



ka 
0± 

yka 



(-9 + 6^/3)Co^l, 



he 



(15) 



B 



ka 



1 + 



-1 + 



r+ 1 



11 13^/3 
16 ^ 32 



(1 _ ^2)5he + 1 ((1 + z^) - (1 - z^){a,a,)t) 



Tihe 
± ) 
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where F^^ is given in (j8| and Zc = 2 — \/3- The two leading terms of near Tc were 
studied before in connection with generahsed extended lattice-lattice scaling [181 [HI [E] • 

2 Scaling theory and CFT predictions 
2.1 Scaling theory 

The singular part of the dimensionless free energjj^ of the two-dimensional Ising model 
satisfies the following scaling Ansatz: 

/sing(5t,5h,K,}) = -9?logbt|-^±(5h/bt|"'^/"%K/l5tr^/''}) (16) 

+ 9hY±{g,,/\gt\y'^/y\{g^j\gt\y^/y^}). 

Here gt, gh, 9u are nonlinear scaling fields associated, respectively, with the thermal field 
r, the magnetic field h and the irrelevant fields The exponent yt is the thermal 

exponent, and takes the value 1 for the two-dimensional Ising model, while yh is the 
magnetic exponent and takes the value 15/8. The irrelevant exponents yj are all negative. 
In the language of conformal field theory (CFT), this scaling Ansatz assumes only a single 
resonance between the identity and the energy. That the dimensions are integers implies 
that there might be multiple resonances which give rise higher powers of log r as observed. 
We have not included such terms in the scaling Ansatz above, as there are no gf{log \ gt\)^ 
terms with n > 1. Following our earlier analysis Caselle et al. |21j discussed the 
scaling theory of the two-dimensional Ising model in considerable depth, in particular the 
conclusions that could be drawn about the irrelevant operators. We discuss this further 
below. 

The nonlinear scaling fields have power series expansions with coefficients which are 
smooth functions of r and the irrelevant variables u = {uj}. In particular one has 

9t = ^a2n(r,w)-/i2", aoiO,u) = 0, (17) 
n>0 

9h = J2^2n+i{r,u)-h^^+\ 

n>0 

n.>0 

In the absence of irrelevant fields, the known zero- field free energy imposes the equalities 
l+(0) = l^-(O) and i^+(0) = y_(0). Furthermore, the known solution for the magnetisation, 

■'in the following, we shall use the notation / — logz — — /3^, with z the partition function per site and 
^ the usual free energy per site. 

''The scaling function y±(a;, {0}), without the effects of irrelevant fields, has been studied recently to 
high precision, see [2Dj and references cited therein. 
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which contains no logarithms, and the known (but not proved) absence of logarithmic terms 
in the divergent part of the susceptibility impose the constraints that the first and second 
derivatives of Y± (0) also vanish. That is to say, 1± (0) = (0) = 0. Aharony and Fisher 
|16] have argued, almost certainly correctly, that there are no logarithms multiplying the 
leading power law divergence of all higher order field derivatives, not just the first two, as 
discussed. In that case it follows that Y± are constants, and further the analyticity on the 
critical isotherm for h requires high- low temperature equality, 1+ = y_. Collecting 
all this information, we have, for the zeroth, first and second field derivatives of the free 
energy, 

f{T,h = 0) = -A{ao{T)flog\ao{T)\+Ao{T), 
M{T<0,h = 0) = Bbi{T)\ao{T)f, (18) 
kBTx±{T,h = 0) = C±{h{T)f\ao{T)\-^-Ea2{T)ao{T)log\ao{T)\+D{T), 

where A, B, C± and E are constants, the background term Aq{t) is a power series in r, 
and the critical exponents are (3 = 1/8 and 7 = 7/4. The free energy and magnetisation 
determine the scaling field coefficients ao(T) and &i(t) which, given our freedom in choice 
of A and B, can be normalized to aQ{T) = r + O(t^) and 6i(r) = 1 + 0(t). The presence 
of any irrelevant scaling fields will manifest themselves as deviations in the predicted form 
of the susceptibility in ( p^ 

To get an explicit expression for the predicted susceptibility in the absence of irrelevant 
fields we start with the zero field magnetization which is known to be the same function 

M = (1 - A;2)l/8 = 2V4A;1/8(i + ^2)l/16(_^)l/8 (^g) 



for all three (square, triangular and honeycomb) lattices. The second equality in (19) 
follows from our temperature definition ([3) and if we use this to solve for 61 (r) in (18) we 
can reduce the zero field susceptibility in (Tsl) to 

fceTxi =C±|r|-^/^F±-Sa2(r)ao(r)log|ao(T)|+Z)(r), (20) 

where 

F± = fci/4(l + T2)V8(^/ao(r))2. (21) 

It only remains to determine ao^r) from the singular part of the zero field free energy for 
each lattice to complete the calculation of F± which we henceforth denote as the Aharony 
and Fisher scaling function F±(A^F). 

It will turn out to be usefuj^ to define the following integral, in terms of which the 
internal energy is defined: 

Iir) = - . — K|— )= , ,x K|^^), 22) 



^According to ([2|, the background contribution D{t) contains terms with arbitrary powers of log|T|, 
which have not yet been interpreted within the context of scahng theory. 

^Identities used can be found in [22], see eqs. 2.597.1, 8.112.3, 8.113.1, 8.113.3, 8.126.3 and 9.131.1. 
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where K is the complete eUiptic integral of the first kind. This function is invariant under 
the high-low temperature change k — t- 1/k. Useful forms at both high and low temperatures 
are obtainable from the Landen transformation, 



K 



2Vk \ 



l + k 



j = il + k)K{k) =(1 + 



K 



(23) 



For the subsequent scaling analysis we will require the singular part of /(r), which is 

2 



It 



sing 



vn/l + T' 
log |r| 



log \t\ • K 



(24) 



2-^1 



2'2' ' 1 + t2 



log |r| • 2-Fi 



1 1 

2' 2' 



We next write the internal energy, per site, in terms of the above integral (22): For the 
square lattice, 

df 

= 2{aiaj)nn = coth(2i^sq)(l - Tl(r)), (25) 

sq 



dK 



where / = —(3^ with ^ the free energy per site. For the triangular lattice. 



5/ 
dK 



tr 



1 + u 
1-u 



3u-l 



2[u3(i_u)3(l + 3u)]i/ 



(26) 



For the honeycomb lattice. 



dK 



he 



J /nn 



1 + z' 
1^ 



l + z 
1 - z 



3/2 



4z - 1 



8[z3(l-z + z2)]l/4 



(27) 



We can calculate the zero- field free energy by integrating these expressions. In fact we 
are only interested in the singular part of the free energy, which we normalise by the 
requirement that it vanishes at T^- With that normalisation, we can write 



dK 



1 1 



(28) 



/sing = - log |r| d-^C,.,F,\^^,-^-l--r 

which is to be compared to /sing = —■^o,q{t)'^ log |r| in (18). The Cj in (28) is the coefficient 
of I{t) in equations (25)-(27) for the appropriate lattice; the dK/dr is also to be evaluated 
with K for the appropriate lattice. 

For the square lattice we determine 2{dKsq/dT)Ci = t/VI + t'^ so that the integrand 
in (28) is seen to be explicitly odd in r and we can identify A^'^ = 1/4 and the even function 



ao(r)' 



Isq 



3 2 , 41 , 

— r H r 

8 192 



147 
1024 



8649 ^ 
81920' 



131072 



(29) 



8 



Equation (29) combined with (21) gives 



1 + -t' 
2 



31 4 125 
384"^ ^ 3072' 



38147 
1474560 



108713 
5898240' 



.10 



+ 0(r 



12\ 



(30) 



which extends the result in [1] to higher order. 

For the triangular lattice we find 8{dKtr/dT)Ci /V27 = r 
^ CnT^""'"^ with the Cn satisfying the three term recursion (n" 

ll/18)c„_i + (n^ -n - 11/18- 15/(144(n2 - n)))c„_2 = 0. The integrand in (g is again 
odd in r and we can identify A^'^ = -v/27/16. For the honeycomb lattice A^'^ 
otherwise the integrand is the same. The scaling functions are 



+ 97rV256 + . . . 
^2 + n + 2/9)c„ + (2n2 



27/32; 



«o(''")^|tr,hc 
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-r H T 
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12- 



and 
F±(A&F)*'"''^^ 



1 + 
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85 



256^ 2048 ■ 
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12^ 



327680 



2621440 



(32) 

We can now compare these scaling functions based on the assumption of no corrections 
to scaling with the observed functions given in (|8|. Define AF± = F± — F±{A^F); then 
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(33) 



The absence of a correction at 0(r'^) in F'^ and Fj^'^ is expected since the operator that 
breaks rotational invariance on the square lattice is not present on these lattices. On the 
other hand Caselle et al. also suggested that because the operator that breaks rotational 
invariance on the triangular lattice first contributes at 0(t®) there might not be any O(r^) 



correction. The clear evidence in (33) of such a correction on the triangular lattice, and 



indeed on all three lattices, shows that there are corrections to scaling operators in the 
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Ising model that are not associated just with the breaking of rotational invariance. We 
elaborate on this in the following section. 

In a similar manner, we can derive the Aharony and Fisher scaling function for the 
kagome lattice. The extra 1 — z"^ term in arises because the magnetization for the 
kagome lattice is given by ([231 eqn. 95]) 



M = Vl - z2 (1 _ 



(34) 



which replaces (19). This introduces an extra \/l — factor into 6i(t) and (21) becomes 

^2 



1 



1 



(35) 



where the denominator in the first term is a normalising factor with Zc = 2 — ^/S the 
critical value on the honeycomb lattice. The remainder of the derivation of the A & F 
scaling function is unchanged, resulting in the relation 



F±(A&F) 



ka 



1 



,F±(A&F) 



he 



(36) 



In view of ( 15 ) and (36 ), we also know that the deviation of the kagome lattice scaling func- 
tion from the corresponding A & F scaling function is identical to that of the honeycomb 
up to a factor, 

pka 1 



AF. 



± 



1 



^AFl". 



(37) 



2.2 Scaling from conformal field theory 

This section draws extensively on the paper by Caselle et al. j21l which was written after 
the appearance of [3]. We adopt the usual notation within CFT. At the critical point, the 
Ising model is describable by the unitary minimal CFT with central charge c = 1/2. The 
spectrum can be divided into three conformal families. They are the identity, spin and 
energy families, commonly denoted [/], [u], and [e] respectively. Each family characterizes 
a different transformation property under the dual and Z2 symmetries. T denotes the 
energy-momentum tensor, so TT is a spin-zero irrelevant operator. Each family contains 
one primary field and a number of secondary fields. The conformal weights of the primary 
fields are hj = 0, /lo- = 1/16, and = 1/2, and all primary fields are relevant. 

The secondary fields are derived from the primary fields by applying the generators 
L-i and L-i of an appropriate Virasoro algebra. L-i plays a particular role, being the 
generator of translations on the lattice, and so gives zero acting on any translationally 
invariant observable. Another important concept is that of a quasi-primary operator. A 
quasi-primary field \Q) is a secondary field satisfying Li\Q) = 0. This condition eliminates 
all secondary fields generated by L_i. As quasi-primary operators are the only ones which 



10 



can appear in translationally invariant quantities, they played a central role in the analysis 
of Caselle et al. j21j . and also in our current analysis, as they are the natural candidates 
for irrelevant operators. 



To make the connection between the scaling Ansatz given in eqn. ( 16 ) and the discussion 



in terms of CFT, we first, for simplicity, set yt to its numerical value, 1, and replace the 



scaling field gt by its leading term r. Then the terms Y± and Y± in eqn. ( 16 ) can be easily 
expanded. They will involve terms of the form 



n(^ =n^ -n^ -ni^r ■ 



As the susceptibility is the second field derivative of the free energy, we must retain terms 
with exactly two factors in the first of the three products above, that is, terms of the form 

|r|2'-i • |r|S''^2 ll\\r\yh \.l\W\yH ' ^ ^ 

II II ig/ M I / jge M I / 

Recall the prefactor ~ before the terms Y± and Y± in eqn. (16). Including this 
prefactor, it is clear that all terms of order in the susceptibility are given by all terms 



in eqn (39) satisfying 

N = 2- (y,, +ya,+J2 p^yi^ + Y.Piy'^^- ^^o) 

The leading term in the susceptibility occurs when there are no e or / fields and yo-i = 2/0-2 = 
yh = 15/8, giving N = 2 — 15/8 — 15/8 = —7/4, which is the well-known susceptibility 
exponent. Exponents for other terms in the table rely on eigenvalue exponents given by 
Caselle et al. [21j, which we summarise in Table [ij 

Caselle et al. [21] have produced a list of irrelevant operators and we reproduce com- 
binations of these operators that contribute to x^'^ ^iid X^^ together with the primary spin 
operator a in Table [2] Power counting as described in |21] and above and leading to re- 
lation ( |40[ ), determines when each combination first contributes. Because corrections at 
O(r^) in both and are not observed and similarly corrections at O(r^) in F^ are 
absent, we adopt the assumption of Caselle et al. that all contributions from combinations 
of the iorn^a'^O^O^iTf)"', n > 0, vanish, as well as all descendants of a'^(TT)'^, and 
consequently these entries are excluded from Table [2] 

Because there are multiple operator combinations at most correction levels in Table 
[2J a unique identification of correction terms with operators is in general not possible. 
Thus the following remarks are to be viewed either as pure speculation or at best a set of 



assumptions consistent with the corrections to scaling that are displayed in ( 33 ) . 
^Here and elsewhere we adopt the convention that is a generic operator in family [x]. 
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Eigenvalue 


Term 


Term 


-2 


QiQi = Tf 


Qi + Qi (sq) 


-3 


Ql + Ql (sq) 




-4 


Qi + Qi (tr) 




-5 


Q| + QI (tr) 




-6 




Qi + Qi (sq) 


-7 
-8 


Q4Q4 


Ql + (sq) 


-10 


Q12 + Q12 




-1 




Qe+Qe M 


1? 


+ (sq) 




-^8 




QIQ'7 + Q7Q3 (sq) 



Table 1: Eigenvalues of various operator combinations that contribute to the susceptibility. 
The spin-zero and spin-12 operators (unlabelled) contribute to both square and triangular 
lattices. The spin-4 and spin-8 operators contribute only to the square lattice (labelled 
(sq)), while the spin-6 operators, labelled (tr), contribute only to the triangular lattice. 



The corrections observed in (33) are consistent with the conjecture that all operator 
combinations of the form O^'O^, O^'O^ or O'^O^O^ are rational multiples of the 
leading order contribution of O"^. Furthermore these multipliers are the same above 
and below Tc- This makes these contributions particularly hard to distinguish from 
the scaling fields associated with the leading contribution. For example, the rational 



coefficient 11/7680 of in AF^^ in (33) is very likely a combination of a direct 
contribution from a'^{Ql + Ql)^ and a scaling field correction from the (t'^{Q\ + 01)^ 
term, whose leading contribution is at order r 



4 



2. We identify all irrational corrections with cr- field operators. Specifically, contributions 
proportional to Cq± with a{Q'^Q'^) and those proportional to Cio± with a{Q'i^Q'^). 



The ambiguity in C%± and Cio± as discussed following eqns. ([8j)-( 12 ) is relevant in the 
present context. A part of Cq± might be a rational number associated with (T'^{Q\Qi) 
and this would further complicate the interpretation of the 11/7680 coefficient of 
described in item[TJ 



3. The coefficients of Cq± in ( |33[ ) on the different lattices are, after dividing out the 
leading term, 1 - 4973t^/5040 + . . . (square), 1 - 403r^/400 + . . . (triangular) and 
1 — 409r^/400 + . . . (honeycomb). Because these are all different, we must conclude 
that the scaling function associated with ^"(QgQs'^) is lattice dependent. An analogy 



is the difference seen in the F± scaling function on the kagome lattice as seen in ( 36 ) . 
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N 


Square 


Triangular 









2 

4 








6 


<y\Q% + Ql? 




8 






10 


crHQi + QDHQl + Ql? 
'^HQi + QifiQiQi) 




12 


a^O^O'' (many terms) 

^(Q| + 0|)'(QfOf) 
{QlQl? 


a^O^O'' (many terms) 

(T{Q?M)iQiQi) 



Table 2: Operator combinations contributing to the susceptibility. The values in the 
first column specify the leading contribution |r|~^/^+^ to x^'i and x*"^ or \t\^ to and 
of the corresponding entries in the second and third columns. 
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It is the equality of F± on the square, triangular and honeycomb lattices to O(r^) 
that is to be considered as "accidental" and not generic. 

4. The very particular structure of the short-distance terms, given in ([2]) is not explicitly 
predicted by CFT. Rather, since the primary logarithm, responsible for the specific 
heat behaviour, is due to a resonance between the thermal and identity operator 
|21j . we might expect additional multiple resonances, giving rise to higher powers 
of logarithms. These are indeed observed, but it does not appear to be possible to 
associate particular operators with these terms — at least not by our naive method of 
just power counting. 

5. Tabled] shows two new distinct cr- field operators at order r^^. If, as we have conjec- 
tured in item[2| each is associated with a new irrational C± then we can no longer 
make any unique identifications as we did for Cq± and Cio±. For all terms in F± 
beyond t^^ we are left only with the numerical coefficients tabulated in appendix A. 



3 Generation of series 

3.1 Quadratic recurrences and Z-invariance 

The algorithm deriving the susceptibility series for the isotropic square lattice Ising model 
with k = sinh^(2/3J), was rather simple using plj^ 

k[C{M, Nf - C{M, N - 1)C{M, N + 1)] 

+ [C*(M, Nf -C*{M - 1, N)C*{M + 1, N)] = 0, 



k[C{M, Nf - C{M - 1, N)C{M + 1, N)] 

+ [C*{M,Nf -C*{M,N -l)C*{M,N + 1)] =0, (41) 

where 

C(M, N) = {ao,oaM,N), C*(M, N) = (ao,o^A/,7v)* (42) 

with the asterisk denoting the corresponding quantities on the dual lattice with the dual 
temperature obtained by replacing k ^ k* = 1/k. Series for the pair correlations can 
be solved iteratively using the series for C(1,0) = C(0, 1) and the diagonal correlations 
C{N, N) and C*{N, N), which in turn follow from the well-known C(0, 0) = 1 and C(l, 1) 
by the Painleve VI type iteration scheme of Jimbo and Miwa [9] , or with a little more work 
from the well-known Toeplitz determinants [28j . 



For the uniform rectangular Isin g m odel, using the methods of their lattice-Painleve III paper 



McCoy and Wu [26] have generalized (41 1 to the so-called A-extended version, in which the coefficient of A-' 



in C{M,N;X) is the j'-particle contribution to the pair correlation function C{M,N). Equations like (41 1 
also exist for n-point correlation functions |27| . 
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Figure 1: Parts of the infinite triangular lattice (black circles), honeycomb lattice (open 
circles) and kagome lattice of rapidity lines (oriented dashed lines) with rapidities u, v and 
w. 

For the isotropic triangular and honeycomb lattices the situation is far more compli- 
cated. We have used the generalization of ( |41[ ) for general planar lattices |24j . together 
with Baxter's Z-invariance \29\ [30] as was first numerically implemented in [31]. More 
specifically, consider the situation in Figure [l| The Ising model on the triangular lattice 
(black circles in the figure) and its dual on the honeycomb lattice (open circles) are Z- 
invariant in the sense of Baxter > with rapidity lines forming a kagome lattice (oriented 
dashed lines) |^ 

To get the isotropic cases we need to choose the three rapidity values as 

u = ^K(fc'), V = \K{k'), w = 0, (43) 

with k' = \/l — k'^ and K{k) the complete elliptic integral of the first kind. 

The interaction constants K = pj are chosen as a function of the two rapidities passing 
through the bond and the directions of these rapidities, following the prescription of Figure 
[2j and as a function of the temperature through the low-temperature elliptic modulus k. 
More precisely, 

sinh {2K{u, v)) = sc{u - v, k') = k-^cs{K{k') -u + v,k'), (44) 
sinh {2K{u, v)) = /c-^cs(u - v, k') = sc{K{k') -u + v,k'), (45) 

where sc{v,k) = sn{v , k) / cn{v , k) = l/cs{v,k). For the dual lattice with k* = 1/k being 

kagome Ising model can be obtained from the honeycomb Ising model by decoration and star-triangle 
transformation [321 IT3] and its spins then live on all the intersections of pairs of rapidity lines. 
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u 



V 



u 



V 



K{u,v) 



K{u,v) 



Figure 2: The two kinds of Ising interactions K and K. On the dual lattice K* and K* 
are assigned similarly, but with modulus k replaced by 1/k. 



the high-temperature elliptic modulus and sinh(2i^*) sinh(2Er) = 1, we have 



sinh (2K*{u, v)^ = k sc{u — v, k') = cs(K(A;') — m + /c') , 
sinh {2K*{u, v)) = cs(n — v, k') = k sc(K(/c') — u + v,k'^ . 



For the triangular lattice we have ( 44 ) with u 



whereas for the honeycomb lattice ( 44 ) with u 



■ V - 

-V 



K{k')/3 or pS)) w ith u - 
■ 2K{k')/3 orl^ withn- 



(46) 
(47) 

2K(A;')/3, 
= K{k')/3. 

Therefore, it is easy to see that the resulting interactions are isotropic for both lattices. 

As the correlation functions only depend on differences of the rapidities, we can add an 
arbitrary common constant to all of them [29]. Changing the direction of a rapidity line is 
equivalent to adding ±K{k') to its rapidity [30]. Together with (43), these two properties 
show that we have invariance under a rotation by 60° for the rapidity lattice, implying the 
required rotation invariance over 60° for the pair correlations on the triangular lattice (or 
over 120° for the honeycomb lattice). In addition we have several reflection properties. 

Most importantly, Baxter's Z-invariance implies that the pair correlation functions, 
apart from their dependence on the modulus k, only depend on the rapidities that pass 
between the two spins [29] , where we have to make all rapidities pass in the same direction 
by adding the above ±K{k') to a rapidity that passes in the opposite direction [SD]. Thus we 
only need to determine universal functions g{ui, ■ ■ ■ , U2m', k) and g*{ui, • • • , U2m', k) giving 
the pair correlations on the lattice {T<Tc) and the dual lattice {T>Tc) These functions 
are invariant under any permutation, or under simultaneous translation by a same amount, 
of all rapidities [29]. As the rapidities Uj can only take the three values (43), we find it 
convenient to introduce the abbreviations 3 1 1 



g[Nu,Ny,Nw] = g{ui,- ■ ■ ,U2m;k) = g[Nw, N^, N^], 
g*[Nu,N^,N^]=g*iui,--- ,U2m;k) = g*[N^,N^,Nu], 



(48) 



■'Compared to [lOj we have interchanged g and g* through this convention. 
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where 

Nu = #{uj\uj = u} , = ^{uj\uj = v}, Ny, = #{uj\uj = w} , (49) 
counting the number of uj^s equal u, v, and w. The symmetry under the interchange of 



Nu and in (48) corresponds to a reflection symmetry that holds in the isotropic case 
(43). Another reflection symmetry gives 

g[M,N,0]=g[N,M,0]=g[0,M,N]=g[0,N,M], 
g[M,0,N] = g[N,0,M], 

g[N,0,0]=g[0,N,0]=g[0,0,N], (50) 

and similar relations hold for g*; these are also reflection symmetries for the uniform 
anisotropic square lattice case represented in Figure [3} 

>< v< >,< >< 

o o ,v o :\ o x o 

• • • • 

o o V o o o 

>' • • • • 

^ o o >' o o o ^ 

>' • • • • 

^ o o >' o \\ o >; o ^ 
uvuvuvuv 

Figure 3: Parts of the infinite square lattice (black circles), dual square lattice (open circles) 
and diagonal lattice of rapidity lines (oriented dashed lines) with rapidities u and v for the 
two directions. 

Now we can invoke the quadratic recurrence relations [23] in the form [10] , 

Sc(n2 — Ml, k')sc{u4 — U3, k') 

x{g{ui,U2,U3,U4, ■ ■ ■ )g{- ■ ■)-g{ui,U2, ■ ■ ■ )g{u3,U4, ■■■)] 
+{g*{ui,us, ■ ■ ■ )g*{u2,U4, ■ • ■)-g*{ui,U4, ■ ■ ■ )g*{u2,U3, •••)}= 0, (51) 

k'^sc{u2 — ui, A;')sc(n4 — M3, k') 

x{g*{ui,U2, U3, n4, • • • )g*{- ■ ■)-g*{ui,U2, ■ ■ ■ )g*{u3,U4, ■■■)} 
+{9{ui,U3,- ■ ■)g{u2,U4,- ■ ■)-g{ui,U4, ■ ■ ■)g{u2,U3,- ■ ■)} = 0, (52) 

where the dots indicate the other rapidities and the modulus that are left unchanged. Eqs. 



(51) and (52) are each other's dual — as is obvious comparing with (44) and (46) — and they 
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can be solved by iteration, once we know the functions g and g* for the two cases with 
ah or all but one of the rapidities equal. Such correlation functions are known under the 
names diagonal and next-to-the-diagonal correlation functions for the square-lattice Ising 
model. An iteration scheme for these is given by Witte ^3], which we adopt with some 

modifications EE 

Let us introduce the abbreviations 
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Xn = (O"0,0fn,n), Vn = (o"0,0O"n,n+l) , = (o"0,0<7n+l,n) , (53) 

for the needed square-lattice correlation functions, together with 

K: = -K{k), £ = -^{k), k={^^''^'^~^' (54) 

where 



Sa = sinh(2Ka) = sc(^|K(A;'),A:'j, Ca = cosh(2i^„) = nci^\K{k'),k'] 
56 = sinh(2K5) = sc(fK(A;'),A:'), C;, = cosh(2Kfe) = nc(f K(A;'), A;' ) , (55) 



for T < Tc, and 

S* = sinh(2Ka) = cs(|K(A:'), A:'), C* = cosh{2Ka) = ns(|K(fc'), A:' 
S; = sinh(2ii:6) = cs(^|K(A;'), A:'), C7* = cosh(2i^f,) = ns(^|K(A;'), A:') , (56) 

for T > Tc- Eqs. (54)-(56) define the rescaled complete elliptic integrals and the hyperbolic 
sines and cosines of twice the horizontal and vertical reduced interaction constants in 
terms of elliptic modulus k. Here Ka is also the reduced interaction energy K^t- of the 
triangular lattice and Kf, is the i^hc of the honeycomb lattice. Duality between low- and 
high-temperature phases is described by the replacements 



k* 



1/k, K* = kIC, £* = k-^ - (1 - A:2)/c) , (57) 

iJb i^a iJb '-'a 

It is easy to check that the dual of the dual gives the original quantities back {X** = X). 



^^In 33^ a single modulus k with < fc < oo is used necessitating definitions such as K< = K(fc) and 
E< = E(A:') for fc < 1 (T > Tc) and K> = K(l/fc) and E> = E(l/fc) for fe > 1 (T < Tc). Here we opt for 
two definitions of k to keep < A: < 1 (cf. (54|). Elsewhere we also adopt a single k; in (3) and (4) our k, 
fcsq, ktr and fchc are to be identified with 1/fcwittc and have been chosen this way because of our numerical 
work for which low-temperature expansions have some advantages. 

^^Here we use the convention of [28] that in (jm,n M is the vertical and A'' the horizontal coordinate. In 
many other works, including fSHI, the opposite convention is used. 
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In general, the nearest-neighbor correlations of the square lattice involve the complete 
elhptic integral of the third kind Ui{n,k) pHllHOllSTl [33] 1 
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yo 



2 Cb 



2 Ca 



for T < Tc, an 



Vo 



Zo 



2_Cl 
vr SI L 

2 c: 



vr S , 



C^Iii{l/Slk)-K{k) 
'C^li^{l/Slk)-K{k) 

C:^Il,{Sf,k)-K{k) 
CfYii{Sl\k)-K{k) 



b ^ 



(59) 
(60) 

(61) 
(62) 



for T > Tc. However, because we have Ka = Ktr, Kb = Khc and the dual/star-triangle 
relation ^ or equivalently 



Ca 

= = CaiCa + So) 



(63) 



these correlations yo and zq (and also i/q and Zq) are also the nearest-neighbour correlations 
of the isotropic triangular and honeycomb lattices. This in turn means they only involve 
the complete elliptic integral of the first kind [Ml EHl [36] . 
To make this more explicit, use [37] 



Hi {-k^sn^{a,k),k) =K{k) 



1 + 



sn(a, k) 



where 



Z(a, k) 



cn(a, A;)dn(a, k) 
Q'ia,k) 



Z(a, k) 



G(a,fe) 



(64) 
(65) 



is Jacobi's Zeta function and [38] 



e{u,k) = e,{z,q)= J2 (-i)V'e 



2inz 



Z 



ITU 



2K{k) 



q = e 



n=—oo 

7rK(fc')/K(fe) 



(66) 



^■^See (4.3a) and (4.3b) of Chapter 8 of [2^, correcting a minor misprint, or (52) of [33], identifying 
ni(n, fe) = n(-n, fc). 

^■*Here, as in (56 1, the asterisk indicates that the RHS is the high-temperature expression. 
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When a is a rational multiple of iK(k'), say a = miK(k')/n, then Z(a, k) can be expanded 
m powers ofgi/'^. It can even be calculated in terms of K(A;), Sa and Sb using the addition 
formula [371 [38] 

Z(n + a,k) = Z{u, k) + Z(a, k) — k^ sn(ii, k) sn(a, /c) sn(u + a, A;), 



and 



TTl 1 TTI 

Z(2iK(A:'),fc) = l{\^{k%k) = -Kl + k) - —. 

Setting u = 2a = 4A,u = a = 2Aovu = a = A = iK(fc') /3 in ([67j), we find 

vri 



Z{A,k) = 
Z{2A,k) 



__ _ ^khn\2A, k) + ^A:2sn2(^, fc) sn(2A, k), 



vri 



Here, using Jacobi's imaginary transformation [38], 

sn{A,k.) = isc(iiK(A;'), A;') = i5a 



sn(2^,A;) = isc(|iK(A;'), A:') = iSb 



1 

kSb 
i 

kSa 



Therefore, 



and 



ICa 

3 

2 a 

3 Sb 



Cb + 



+ 



+ 



Cb _|_ 1 Cg ^ 1 SbCg 

Sb 2 SaSb 6 S^ 

Ca 1 Cb 



•^a 3 5",^ 



/C, 



Ca ^ 1 Cb 1 5bCf, 



SaSb 2 Sb 6 5"^ 
Cf, 1 SbCa 



SaSb 3 (S,^ 



/c, 



for T < Tc, 
/C, 

for T > Tr. 



(67) 
(68) 



(69) 

(70) 
(71) 



(72) 



(73) 



Results (72) and (73) differ by duality as defined in (57) and (58). 

We next rewrite ( 72 ) and ( 73 ) using ( 63 ) or alternatively using Ka = i^tr and Kb = i^hc 
with the explicit connections to u and z given in (|4| and ([5])(^ With the latter we obtain 



2/0 



l + u 
3(1 -u) 

l + z^ 
3(1 -z2) 



1 + 



2 + 



2(1 - 3u) 



/C 



V(l-n)3(l + 3n) 
(l + z)(l -4z + z2) 



{l-zf 



IC 



(74) 
(75) 



Cf. also (1061 and (1071 in the following section. 
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which can be compared directly with the internal energy results in Table I of Houtappel 
1 34 1 PI These results are also the basis for our p,6n and i27n- the equality follows by using 
(22) and the low-temperature Landen transformation from (23) to yield /(r) = 2^/klC. 

We can also rewrite (44) and (54) of [33J. Then the first few square-lattice correlations 
in the low-temperature phase are 





1, 




Ca ( 


yo = 


35,1 




Cb 1 


zo = 


35fe \ 


yi = 




Zl = 





Xl 



l-{Ca-2Sa){Ca + SaflC], 
{Cb-2){Cb + l? 



2 + 



SI 



)yo + , 

>->a ' 'Ja 



whereas the corresponding quantities in the high-temperature phase are 

Xq 

Vo 



^0 

y*i 



a 1 



Xl — C ^ 

iCb-2){Cb + lf 



JC 



Cj2 + iCa-2Sa){Ca + SafjC* 



Sb^ 
Sa 

Sb 



SbCa 



£8* 



(76) 
(77) 

(78) 

(79) 

(80) 

(81) 
(82) 

(83) 
(84) 
(85) 



These results are fully consistent with duality defined in (57) and (58). In addition, we 
have Zq = Ca — Sayo and Uq = Cb — SbZQ, in agreement with (11) in [21]. 
Witte's initial conditions, (40) and (42) in [33], can be rewritten as 



ri 



ro = 
2k £* 



n 



£* 
~£ 



(86) 
(87) 



^The results of Wannier [35] and Newell [36] differ by Landen transformations |38| 



2Vk 
1 + fc' 



■•^Wannicr 



l-k' 

l + k'' 
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and 



"•; = 1, >=; = 1, 

Then further quantities can be found systematically using 

(2j + 3)(1 - rjr,)rj+i = 2j(k + k'' + (2j - l)r,f^-i^ 

-(2j - 3) (l + (2j - l)r,f,)r,_i, (90) 

(2j + 1)(1 - r,f^)f,+i = 2j(k + fc-i - (2j - 3)f,r,_i)f, 

-(2j-l)(l-(2j + l)r,f,)f,_i, (91) 



and the identical equations for r* and r*, see (38) and (39) in [33]. The further diagonal 
and next-to-the-diagonal correlations follow using 



-1 










(1- 












Xj+l 






Sb 


Xj 








Xj + i 




rj+i 


Sa 


Xj 






Sb 



(92) 



and their dual versions obtained by replacing all quantities by their * versions. These 
last few equations can be found combining (31), (36), (59), (63) and (64) of [33j. For the 
current purpose one only needs and z* for n = 0. 

We have now all equations from the square-lattice Ising model needed to generate g and 
g* with all or all but one of the rapidities equal in a form that makes the lattice symmetries 
and duality manifest. Thus we can now construct a "polynomial-time" algorithm for the 
high- and low-temperature series coefficients for the susceptibility of the isotropic Ising 
model on triangular, honeycomb (and kagome) lattices. For efficiency of the algorithm, we 
desire series with only integer coefficients. Series in the low-temperature u = exp(— 4i('tr) 
are certainly acceptable; because the coefficients in these series can be reduced to lattice 
counts, they are necessarily integer. A useful alternative in the square lattice case [4J 
was the elliptic parameter k. The corresponding alternative here suggested by the kt^{u) 
relation Q is an expansion in the variable k = (A;^/16)^/'^. Inversion of kij-{u) results in 
the series 

u = k-2k^ + h* + 3P - IQk^ + + 40fc8 - I61fc9 + l^fcio + . . . (95) 

3 9 81 
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and although the rationals in (95 ) can be ehminated by the change k ^ k/3 the coefficients 
in any correlation function series in k will grow unacceptably rapidly. A third alternative 
is expansion in g^/^ where q is the elliptic nome. This is suggested by A: = (A;^/16)^/^ and 
the known expansion k^/16 = q-Sq^ + .... 

All our elliptic functions have natural expansions in terms of the elliptic nome 



using Jacobi theta functions, i.e. 



i(0,g) 



^3(0, 



exp 



k' 



7rK(A:') 



1(0, g) 
3(0, (?) 



/C= [03(0,(7)]', £=W,q)f 



t(0,g) MO,q)Y 



Sa = -i Snf ii K{k'), k), Ca = CYi(\\K{k'), k) , 



Also, from (55), 



Si, = -i sn(^|i K{k'), kj, Cb = cn(^|iK(A;'), kj , 
using Jacobi's imaginary transformation. In terms of theta functions, 

-i Oi{Za,b,q) ^ _ . d2{Za,b,q) 



Sa,b 



with 



Vk 0A{Za,b,q) 

vr iK(A;') 



Ca,b 



k 9i{za,b-,qy 



Zb — 2Zq,, 



vV6 



" 2K{k) 3 

From the above we expect to end up with expansions in the nome 

7rK{k' 



exp 



3K(fc) 



1/3 



(96) 

(97) 
(98) 



(99) 

(100) 
(101) 

(102) 



and this is the good expansion variable that we 

usedj^ 

For expansions in terms of the nome it is also advantageous to break the symmetry 
defining 

rj = i-kVPj, fj = i-kr^p,, (103) 
and similar for and fj, in order to avoid square roots of the nome. 



^^There are many other cases where series in the nome are advantageous. Whenever ah rapidity diflterences 
are of the form mK{k')/n with fixed integer n, we can expand the susceptibihty in powers of g = q^^", see 
the text following 1 66 1. 
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3.2 Alternative expressions 

The functions sc(|K(/c'), k') and sc{^K{k'), k') have an algebraic representation in k which 
one can obtain by expanding identities such as cs(gK(/c) + |K(/c) + ^K{k),k^ = using 
standard addition formulae and then solving the resulting quartic equation for sc(^K(/c), k). 
One finds 

where 



sc(^K(A;'),A:' 



(104) 



R 



X + y/3-X^ + {k-^ + k)/X, X = ^l+((A:-i-A:)2/4)'^'. 



(105) 



Note that R is self-dual, i.e. invariant under the replacement k ^ 1/k, while sc(^K(/c'), k') 
■H- cs(|K(A;'), k'). If we take k = ktr, the low-temperature elliptic parameter Q, then one 
can verify 



1 



1-u 



= sinh(2i^i 



tr) 



and 



rR _ ^(l-u){l + 3u) _l-z^ 



sinh(2Khc), 



(106) 



(107) 



k 2u 2z 

where in (107) we have used ([s]) for u{z). In this way we confirm directly from (104)-( [T07 ) 
and the definitions (55) that Ka = i^tr and Ki, = K^^- 



The expansions in the (cube root) nome q = exp(— vrK'/SK) described in the preceding 
section can be applied to (C^ — Sa)"^ to give directly u = u{q). We obtain the formula 



u 



oo 



4n 



7 8n+2\ 



n=0 ^ ' ^n=0 

2q^ + 3^^ - 4f + 7q^ - 12g^^ + 17^^^ 



oo \ 4 

3n.(n+l) 



24gi5 + ... 



(108) 



which explicitly shows u{q) is an integer series. Whether correlation function series in u or 
q will show the slowest growth in the magnitude of the series coefficients depends on the 
singularity structure of the correlation functions. Now the correlation functions as series 
in u have radius of convergence 1/3 governed both by the ferromagnetic singularity at 
M = 1/3 and an unphysical singularity at u = —1/3. There are other more distant complex 



singularities and what we have found numerically and describe in section 4.2 is that there 



is a close analogy with the singularities on the square lattice. Indeed we conjecture that 
\ktr\ = 1 is dense with singularities and parl^^ of a natural boundary for the triangular 
lattice. Now the circle |A:tr| = 1 maps to arcs in the g-plane with distance to the origin 
bounded below by exp(— 7r/3) = 0.3509. . . and it is this distance that fixes the radius of 
convergence of the q series. It implies that asymptotically in A'^ we have terms of magnitude 



'For the complete natural boundary see Figure [s] in Section 4.2 
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~ 2.85 q compared to As an example of what we observe in practice, the 

coefficients in the series expansion of the low-temperature triangular lattice susceptibility 
are, at the largest N we have available, dominated by a single u-plane singularity pair 
(u — {1 i: 2i)/5) giving a q series coefficient dependence of mag nitude ~ 2.78^ /N^^/"^ . 
In conclusion, there is coefficient size reduction in going from series in to g but it is not 
dramatic. 



3.3 Computational details 



As discussed in subsection 3.1 the calculation of the triangular and honeycomb lattice 
susceptibilities as high- and low-temperature series of length requires as an intermediate 
step the calculation of two triply indexed arrays g and g* . That is, we are dealing with 
0(A^^) elements, each element being a series of length A^ with integer coefficients whose 
(digit) size increases linearly with A^. Fortunately this 0(A^^) memory requirement can be 
circumvented by a careful sequential arrangement of the calculation and the description 
of this with emphasis on the storage structure we have implemented is the content of this 
section. 

To begin the discussion we show in Figure |4] the i,k,j triples indexing the correlation 
functions C{R) = g{i, k,j) at the triangular and honeycomb lattice sites on the minimum 
sector necessary for obtaining the susceptibility. To obtain these triples refer to Figure 
[l] and let be the number of rapid ity lines of type a between the origin and site R. 



Then the rules given in Section 3.1 can be summarized by saying that if R is in the 
7r/3 sector above(below) the horizontal through the origin then g{i,k,j) = g{Ny,Nu,Nu]) 
(= g{Ny, Nyj, Nu)). Clearly if R is on the horizontal, = and A„ = A^. If at 
least one of i,k,j is zero then the corresponding g is also a correlation function on the 
anisotropic square lattice. For example, if i = 0, then in g(0,k,j) we can identify j = Ux 
and k = Uy where and Hy are the Cartesian coordinates of sites in the ffist quadrant 
of the square lattice obtained by rotating that shown in Figure |4] clockwise by 7r/4. In 
the fourth quadrant of this rotated lattice g{i,0,j) = g{—ny,0,nx)- With the exception 
of g(l,0, 1) only elements g{0,ny,nx) are needed as initial conditions for the recursion 
relations for the general C{R) on the triangular and honeycomb lattices. 

It is worth remarking that the correlations g{0,n,n) on the triangular lattice diagonal 
can be calculated as Toeplitz determinants [39] and we have used this as an important 



check of our computations. We also note that of the symmetries (48) and (50) satisfied by 
the g{i,k,j), two in particular that we use below are g{i,k,j) = g{j,k,i) and g{0,k,j) = 
9(0, j,k). 

One observes in Figure |4] that within each "shell" Ng, that is, sites between lines A^s — 1 
and A^s, the central k index is either 2A^s — 2 or 2A^s — 1. It turns out that with a few 
exceptions, an array at fixed k can be computed from elements in an array with index 
k — 1. Proceeding sequentially through "shells", or equivalently k, reduces the memory 
requirement from 0(A^^) to 0(A^^). It also has the advantage of allowing the calculation 
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Ns=2 / / 



Ns=3 




Figure 4: A vr/S section of the triangular/honeycomb lattice. Small triangles mark the 
triangular and even honeycomb lattice sites. Numbers are the i,k,j indices of the g and 
g* correlation functions. Lines labeled by Ns are the boundaries of the "shells" described 
in the text. The inset shows the g and g* labeling of sites on the anisotropic square lattice. 
Vectors are the directions of the rapidity lines on the two lattices. Dashed lines labeled x, 
y and z indicate the diagonal and near diagonal elements calculated by the Witte recursion 



relations described in section 3.1 Note that the two lowest rows on the square correspond 
exactly to two rows on the triangular lattice. 
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to be stopped and restarted if necessary at convenient intervals and makes calculation with 
an N of several hundred to a thousand practical. 

We take the arra}^^ g(i,k,j) = gk{i,j) indexed — using Maple notation — in the double 
sequence seq(seq(gk(i,i+2*j) , j=0. . .Ns-i),i=0. . .Ns) for even k = 2Ns — 2 as con- 
stituting "shell" Ns{a). The array g{i,k,j) = gk{hj) with odd k = 2Ns — 1 and indexed 
as seq(seq(gk(i,i+2*j + l) , j=0. . .Ns-i) ,i=0. . .Ns) constitutes "shell" Ns{b). The in- 
dexing for both arrays is such that j > i,i + j < k + 2, and satisfies the requirement that 
i + k + j be even. As a specific example of this indexing, the A'^, = 3 case is illustrated as 



(109) 



(54) 
(55) 
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(3a) 

{3b) (109) 



with gk label on the left and "shell" label on the right. Both Ns{a) and Ns{b) arrays are of 
length L = {Ns + l){Ns + 2)/2 and we introduce a third notation, namely the single indexed 
gki^),(^ = 1 ■ ■ ■ L. The physical elements on the lattice are only a subset of 3Ns — 1 elements 
in each array. Specifically, the triangular and even honeycomb sites are at locations £ = 
L — 1 — n{n + l)/2, n = 1 . . . Ns and are indicated by the double arrows in (109) 
honeycomb sites below the horizontal in Figure |4] are at £ = L — n{n + 1)/2, n 
while those above are &t £ = L — 2 — n[n + l)/2, n = 2 . . . Ng. Both sets are indicated by 



The odd 
... A^,, - 1 



single arrows in (109). 



We also require linear arrays which for identification purposes we will denote as with 
the even and odd k arrays being distinct. The array do is initialized by elements from the 
anisotropic square lattice array x described in Section 3.1 , di by the corresponding elements 



from y. In subsequent calculations, dfc-2 will be renamed dk and certain elements changed 
by an in-place replacement determined by the quadratic recursion formulae. Details will 
be described below; for now it is enough to know that the changes will maintain dk{l) = 
g{0,k + 2,k mod 2), 4(2) = g{0,k,k + 2) and 4(3) = g{0,k,k + A). The4(n),n > A^^ + 1 
remain unchanged from the initializations 



(do) 020 002 004 006 008 



(4) 031 013 015 017 019 



iao) 

{91) 
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fr 
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(16) (110) 



where arrows indicate physical site elements as in (109). The underlines indicate elements 
that have been copied, specifically do{n) = Xn-i and di(n) = Un-i for n > 1. Also, the 



While we only refer to an array g here, there is a strictly parallel dual array g* . It is to be understood 
that such convention applies throughout this section. 
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elements in go are xo, xi and zq while the first two in gi are yo and yi. A special remark 
is in order for elements do{l) and di{l) — these are equal respectively to do(2) and di{2) 
because of the symmetry gk{0,j) = gj{0, k). The third element in gi is given by 



g{l, 1, 2) = g{0, 1, 1)5(1, 0, 1) + /(0, 1, 1)(5*(1, 0, 1) - g*iO, 0, 2))k- 



tr 



(111) 



which is a special case of the recursion equation ( 116 ). Note that the dual of ( 111 ) requires 
both g ^ g* and /ctr — 1/A;tr- 

This completes the initialization except for combining the physical site elements in 
(110), with appropriate multiplicity factors, into (summed) correlation functions from 
which susceptibilities will be determined as a very last step. These functions are chosen to 
distinguish between even and odd sites; given the initialization (110) we set 

Ce = 5<7o(l) + 6<55i(l), 
C:=gl{l) + Q5g\{l), 



Co=355o(3) + 6(^<7i(3), 
Q*=35S(3)+65l(3). 



(112) 



Each 6g in (112) is the magnetization subtracted g — M which applies only to the low- 
temperature variables and not the high-temperature duals. 

The recursion in which new gt are calculated starts with k = 2 and Ng = 2. In 
the general case the first element of g^ is initialized by copying from (ifc_2, specifically 
(7fc(l) = dk-2{^) = 9(0, k,k mod 2). We then proceed sequentially from the gki2) to 
the final gk{L),L = {Ng + l)(-^s + 2)/2, using the quadratic recursion relations for each. 
Unless forced otherwise, we use only elements from gk and gk-i to minimize what is kept in 
memory and this requires that different forms of the recursion relations be used depending 



on the i,j combination in gk{i,j). In the order used, these are^*^ 

- [5fc-i(0,i-l)' + (5Li(0,i-l)' 

-gl{0,j - 2)gl_^{0,j))Rhr]/gk-2{0J - 2 
: [gk-i{0,k + lf + {gl_^{0,k + lf 
-gl{0,k)dl_2i3))Rhr]/gk-2{0,k), 
: [5fc-i(0,l)2-(5Li(0,l)' 

-gU0,0)gl.2ihl))Rktr]/9k-2i0,0) 
■■ [5fc(0, j - l)c/fc_i(l,i - 1) 

HgWJ - l)9*k-i{hj - 1) - 9m, J - 2)gl_,{0,j))ktr 
/5fc-i(0,i-2), j>l 

■■ [5fc(i - l,i - l)5'fc-i(« - l,j) 



gk{o,j) 

gk{0,k + 2) 
9k{l,l) 

gk{i,j) 



l<j<k, (113) 
(114) 
(115) 



gkihj) 



+i9kii - 1. J - - 1. i) - 9*k{i - 2, j)fffc_i(i, J 

/5fc-i(i-2,j-l), i>2 



(116) 



(117) 



As noted in the context of (|109| each equation is to be understood as a pair. Here the second member 
is obtained by the interchange gl and replacements d* ^ d and fctr — > 1/fctr- 
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where the (self-dual) multiplier R is given in (105) or, more simply, as R = Sb/Sa by 



combining (106) and (107). The j index in these recursion equations increments in steps 



of two to maintain i + j + k even, a condition that also eliminates (115) unless k is even 



Indexing functions are easily established which relate the location of the right hand side 



elements in (113)-(117) to those on the left; this is a coding detail that we do not give 

gk{j^i) may have to be invoked to 
2, A; + 2) which in 



here except to remark that the symmetry gkihj] 
locate an element. The special element d^_2(3) in (114) is g*{0,k 



our construction of the gk-2 array was explicitly excluded from being one of the elements. 
As an observation on memory requirements, only the first A^^ elements of array gk-2 are 
required for implementing (113)-(115) so that most of the memory used by gk-2 could be 



released before the gk calculation is started. For all further calculations in (116) and (117) 
only gk-i need be maintained in memory. In fact with a small location offset of 2Ns + 1 the 
replacement gk-i — ^ gk could be done in-place and thus reduce memory requirements even 



further. On completion of the gk calculation in (113)-(117) the gk elements corresponding 



to physical lattice sites are accumulated into the C and C* as in (112) with appropriate 
attention to multiplicity. 



We must also update the dk-2 array that has just been used in (114) in preparation 
for su bsequent iterations in k. The dk, and for completeness the relevant gk, are shown in 
(fTl8) 



(do) 
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(ffi) Oil 013 



t 

101 (la) 

t 

n? (16) 



(d2) 040 024 026 006 



(32) 020 022 024 121 123 222 (2o) 



(d4 



(53) 031 033 035 132 134 233 (26) 
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fr 



051 035 037 017 

t fr t 
(34) 040 042 044 046 141 143 145 242 244 
060 046 048 028 008 ... 

t t t 



343 (3a) 



it 



3) 051 053 055 057 152 154 156 253 255 354 (36) 



071 057 059 039 019 

t t 



(ge) 060 062 064 066 068 161 163 165 167 262 264 266 363 365 464 (4a) 

(118) 

to illustrate the changes in the dk as one proceeds through to the completion of "shell" 



3 and into "shell" 4. The notation in (118) is as in (109) and (110). Of special note are 
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the underlined elements in dk,k > 2, and the fact that all changes are made in-place. 
Specifically this means that we first rename dk-2 to dj.- Then we copy dk{Ns + 1) to 
dfc(l) since it is needed both in its original location where it will be overwritten and in 



a subsequent gk+2 calculation The second copy is from the just completed gk{Ns + 1) 
to dk{2). The transformation of dk is then completed by a sequence of in-place quadratic 
recursion transformations of elements dk{n) starting at n = A'^, + 1 and decrementing to 
n = 3. Each recursion is given b} 
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dk{n) = [dk-i{nf + (4_i(n)2 - ^(n - l)4(n + l))RK]/dk{n) (119) 
which is in a form identical to ( 113| ) including R from ( 105 )-( 107). All operations in "loop" 



k have now been completed and we can restart the overall cycle begun following (112) after 
incrementing — )• A; + 1 and, if the new k is even, Ng ^ Ng + 1- 

On completion of all recursions, high- and low-temperature susceptibility series are 
generated from the C and C* as follows. The triangular lattice susceptibility for T < Tc is 
given directly as 

k^T x'l{u) = C,{u) (120) 

while that for the honeycomb follows from the duality /star-triangle transformation ([s]) and 
is 

k^T x^-{z) = Ce (u = z/{l -z + z^)) ±Co{u = z/{l -z + z^)) . (121) 



Note that both odd and even sites contribute in ( 121[ ) with the sum for the ferromagnet; 



the difference for the antiferromagnet. The results for T > follow by duality and are 

k^TxXiv) = Cl{u = v/{l-v + v^)), (122) 
k^Tx'^^iv) = Cl{u = v^)+C:{u = v^), (123) 

where v = tanh(/C) is the conventional high-temperature variable and K is Ktr or ii'hc 
as appropriate. All of these susceptibilities agree with the earlier work by Sykes et al. 
[IQIIIIIE]. 

In our implementation of the above procedure we made full use of Maple's automatic 
series multiplication routines in full integer arithmetic. This is similar to what was done in 
[1] for the square lattice and allowed us to reach series of adequate length. However we did 
introduce several modifications to improve efficiency. First, as also in f3], when generating 
high- and low-temperature series the recursions were set up to deal directly with the much 



^^There are in general other elements that could be saved for gk+2m,'m > 1, but we have opted instead 
for a small amount of redundancy in our calculation. 

^■^Once again there is a second member obtained by d -o- d* and kti — )■ 1/fctr which must be done before 
n is decremented. 
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smaller residuals 6g = g — Af^. As an example of this change, the recursion (117) becomes 



^dkiij) = Sgk{i-l,j -l) + 6gk-i{i-l,j) -Sgk-i{i-2J -I) 
[{6gk{i - l,j - 1) - 6gk-i{i " 2,i - 1)) 
x{Sgk-iii - 1, j) - Sgk~i{i - 2, j - 1)) 

+ {gl(.i - 1,J - l)9l-i{i - - glii - 2,j)gl_,{i,j - l))kt,] 
/{M^ + 5gk-iii-2,j-l)), i>2 (124) 

in which the magnetization appears only in a denominator factor. 

A second change was based on the observation that all g* terms on odd honeycomb 
sites are of the form ^/u times series in u. If we define these g* terms as g*ktr and use 
g* in the recursion relations in place of g* one can eliminate all occurrences of ^/u and 
dramatically speed up Maple's handling of the resulting series 

Thirdly, we transformed from series in u to series in the (cube root) nome q = e~'^^'/3K^ 



As remarked in Section 3.2, the effect is not dramatic but because the implementation of 
a variable change is so easy we did take this opportunity for improved efficiency. 

For the high- and low-temperature susceptibilities, we generated series to "shell" 160 in 
about 40 days on a 3 Ghz Pentium processor with 500 Mbyte memory. This gives x^'^iu), 
X^'^iv) and x^'^i^) to about 640 terms and x^^i''^) to about 320 and these series can be found 
in p]. 

We have also run the recursion program for series in r to 0(r^^) for the data necessary 
to determine the "short-distance" terms in x- Here there is no magnetization subtraction; 
instead 5*(t) = g{—T) and the code simplifies considerably. It is only practical to run in 
floating point and we have gone as high as 121 "shells" with an accuracy estimated better 
than about 500 digits. Another difference from the high- and low-temperature series case is 
that the correlation data from different shells is not accumulated but rather kept separate 



so as to allow a fitting procedure completely analogous to that described in [J, Section 
6]. Note that there is a distinction between what constitutes a shell for short-distance 
fitting and the "shells" as defined in Figure |4j First, a fitting shell contains only one layer 
each of odd and even sites — not the two shown in Figure |4] — so our data extends to 241 
fitting shells. Secondly, we try to keep fitting shells as close to perfect hexagons as possible. 
Symmetry dictates that we use even sites from a single gk but odd sites are taken from gk 
if they are below the horizontal in Figure |4] and from gk+i if they are above the horizontal. 

The individual shell values are of little intrinsic interest and are not recorded here. 
Instead we give "short-distance" terms, which are the output of the fitting procedure, in 
an abbreviated form in Appendix C and to the full estimated accuracy of our calculations 
in |43j . To complement the much longer 2042 term high- and low-temperature square 

^^The rescaling also has the advantage of eliminating about one-half of all explicit occurrences of fctr in 



the recursions (113l-(117l and thus reducing the number of required series multiplications 

^''One important observation is that the factor ^/s that appears in various equations in [4] is now to be 



interpreted as fc^^* — it remains as the same function of r. 
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lattice series from |44] . we have rerun the code in [4j for series in r to 0(t^^) to 241 shehs. 
Our extended fits confirm the earher results from [4J and the new output is recorded in 
Appendix C and [43] as for the triangular/honeycomb data. 



4 Extracting the scaling function 

4.1 Changing the series variable 

Once we obtained the high- and low-temperature susceptibility series, we analysed them 
to extract the scaling function. Firstly, we normalised the series variable so that the 
ferromagnetic singularity occurs at 1. For example, for the high- and low-temperature 
square lattice series we use the variables z = s and z = respectively. 

We began with short-distance terms calculated from the expansion of the susceptibility 
in terms of r as described in [U Section 6], and a number of Aharony and Fisher scaling 
terms which are known to be accurate. These are of the form T"(ln|r|)^ and r"'^/^"'"" 
respectively. 

We converted these to series in our chosen variable z in the following manner. First 
we expressed each of these terms as a series in 1 — z, of order approximately 50, which 
may be multiplied by (ln(l — z))'' or (1 — z)~'^^^. Each term in the 1 — z series was then 
expanded as a series in z to the full length of the susceptibility series (about 2000 terms for 
the square lattice, for example), and the results added up to produce a series in z for each 
short- distance and Aharony and Fisher scaling term. All these series were then subtracted 
from the susceptibility series. This formed a new series J2n (^nz"", singular at z = 1. 



4.2 Singularity suppression 

The next step involved suppressing the effect of the competing singularities on the series. 
For the square lattice, the singularities of the susceptibility are given ([Sj) by the singu- 
larities of the A^-particle contributions. These lie on the unit circle |s| = 1 at the points 
Ski = exp(20), where 

2TTk 27r/ 

2 cos (9 = cos— -h cos— , < /c,Z < iV (fc,/ not both 0). (125) 

For the low-temperature series, only the even-A^ singularities are relevant. The asymp- 
totic behaviour of the susceptibility near each of these singularities is given (jUj) by 
(1 — z/z'Y, where z' is the singularity and p = (N"^ — 3)/2. This introduced a term 
into the susceptibility series which behaves asymptotically as n^'^^^ (since \z'\ = 1). This 
has the potential to dominate the effect of the scaling term 7-"^/^+'^ ~ (1 — 2;)"'^/^+'^ for 
large a, since it introduces a term into the susceptibility which behaves asymptotically as 

j^3/4-a 
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The simplest procedure (which was the one used in [4j) to rectify this is simply to 
multiply the series by 1 — zj z' . This changes the behaviour of the contribution from the 
singularity at z' to n~P~^, but leaves the contribution from the scaling term at n^/^"*^. 

However, because we know the exact form of the singularity, we can use a more accurate 
suppression. To illustrate, we begin by observing that 



1 



+ 



p + 1 



dz 



1 



c, 



(126) 



where c is a constant. Expressing the original singularity term as a series shows that 
applying Ip forms the new series 
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P + 1 



n 



(127) 



which completely removes the (1 — z/ z'Y term. Moreover, because 



Z \ P+a 



{p + a + l) 



Z \ p+a+l 



+ C, 



(128) 



this transformation also has the additional effect of suppressing (1 — zjz'y^^. In other 
words, the contribution to the susceptibility from this singularity goes from n^^^^ to n 



-p-3 



when we apply this suppression, compared with n ^ ^ when we simply multiply by 1 — z/z'. 
In addition, applying the integral operator to scaling terms gives 



/,(l-z)-V^+^=(l-^) [l-zY 



7/4+a 



P+l 



z'(-7/4 + a + l) 



:i 



-3/4+a 



(129) 



which still contributes Tn?/^~'^ to the asymptotic behaviour of the susceptibility series. So 
this operator suppresses the competing singularity while not asymptotically affecting the 
scaling term. 

An unfortunate consequence of applying Ip for a complex singularity is that the series 
resulting from (127) has complex coefficients. This can be avoided by observing that since 
the susceptibility is real, for every singularity z' there is a corresponding singularity of 
the same order at z' . Sequentially applying the suppression to both of these singularities 
results in the series 
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(130) 



which can be seen to have real coefficients. As all we are doing is applying formula (127) 
twice for two different singularities, the effects on the singularity and scaling terms that 
we observed above still hold. 
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In practice, we also suppress the higher-order terms (1 — zj z'Y^'^ for o = 2, 4, . . ., using 
the above suppression formula (with p replaced by p + a) for each a. The maximum a that 
we use varies for each singularity and is determined empirically as described below. 

For high-temperature series, only the odd-A^ singularities are relevant. The asymp- 
totic behaviour of the susceptibility near each of these singularities is given ([5]) by (1 — 
zj z'y ln(l — zj z'\ where "p = (A^^ — 3)/2. These terms can also be suppressed by the same 



formula (127). This can be seen to be true because applying the same integral operator 
results in an analytic term for integer p (which is true for odd A^). Again, we suppress a 
number of higher powers. 

In order to determine which singularities should be suppressed and by how much, we 
apply a Fast Fourier Transform diagnostic, as described in [44 ^ Section 7]. We first do 
a preliminary fit of the series to our functions, as described in section |4.3| below, and 
subtract the fit from the series. The dominant unsuppressed singularity in the remainder 
is expressed by periodic behaviour of period 27r/^, for a singularity located at exp(z0). By 
applying FFT to the remainder, we can observe the periods of the dominant unsuppressed 
singularities, match these to the known singularities, and increase the suppression on these 
singularities (by suppressing more higher-order terms) . This is repeated until the remainder 
has a satisfactorily small amplitude. 

The analysis of the triangular and honeycomb series is almost identical, though we 
must suppress the appropriate singularities (see 06]). For these lattices, it is conjectured 
that the singularities lie on the curve of Matveev and Shrock ( [U] ) , 

3 

1 + 3m^ - 2u(l - u)x = 0, --<x<3. (131) 

We further conjecture that the singularities are given implicitly by this equation when x 
takes the values 

27rfc 27r/ 27rm 

Xkim = cos + cos + cos ^ , + £ + m = mod A/ {k,l,m not all 0). (132) 



To suppress these singularities, we again apply formula (127), assuming that the form of 
the singularities, and in particular the exponent p = {N'^ — 3)/2, is the same for these 
lattices as for the square lattice. 

Partial confirmation of this conjecture arises from the singularities that we observe 
from our FFT diagnostic as we suppress singularities. We have observed the singularities 
corresponding to this formula for (k, I, m) = (1, 0, —1) for A^ = 3 to 8 and A^ = 10, and for 
{k, I, m) = (2, 0, -2), (2, -1, -1) for A^ = 6. 

We checked for additional singularities by analyzing both low- and high-temperature 
series of the triangular lattice susceptibility in the (cube root) nome q = exp{—TrK'/3K). 



Because all complex portions of the g-plane curves defined by (131) are at or within the 
distance to the ferromagnetic singularity, the high order series coefficients in q will be 
dominated by the complex singularities. By a succession of suppressions of the dominant 
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terms and FFT diagnostics we have identified the same = 4, 6 and 8 singularities as 
found in the ti-plane; in addition {k,l,m) = (2,-1,-1) for = 8 and two singularities 



consistent with x ~ —1.21 and -1.35 in (131). The latter singularities are those on the 



left, upper plane arc shown in Figure [5} From the T > Tc series in q we find the same 
= 3, 5 and 7 as in the t;-plane, the {k, I, m) = (2, —1, —1) for both N = 5 and 7 and a 
singularity consistent with x ~ —1.27 in ( |131 ) and shown on the left, lower arc in Fig ure [S] 



The singularities on the left arcs are not identifiable with any small integer values in (132) 



Thus, although we propose that the closed curve in Figure [5] is a natural boundary for both 



low and high temperature, we can only give (132) as the conjectured singularities for the 



right arcs corresponding to x > — 1 and |A;tr| = 1- Confirmation of this and a formula for 
the singularities on the left arcs can presumably be obtained by an analysis of the Vaidya 
(|46j) integrals. 

4.3 Fitting 

Once all the singularities are suppressed, we fit the series to our scaling functions. We 
use only terms which are known (or assumed) to be nonzero, and leave out the known 
(and removed) Aharony and Fisher scaling terms. In other words, we fit to the linear 
combination 

iVlT^ + t)1/V-^/4 (aer^ + agr^ + aior'' + . . .) (133) 

with og, as, aio, ■ ■ ■ our fitting coefficients. 

Firstly, we convert each term in this expression from r to our series variable z, as 



described in section |4.1[ We then apply the singularity suppression that we applied to our 
susceptibility series to these fitting functions, so that the required equality between the two 
functions is maintained even though both functions have been changed by the suppression. 
Finally we fit the suppressed series to the linear combination of our suppressed fitting 



functions. Suppose that the transformed and suppressed fitting function (133) is Y2n fnZ^ , 
while the subtracted and suppressed susceptibility series is X]n'^«^"- W^e choose the am- 
plitudes to minimise the expression 

"■2 

Y^^fr^-Cnf. (134) 
n=ni 

The range of n in the sum can be varied, but we always choose n2 to be the largest 
available power of z in our susceptibility series. In addition, varying ni will change the 
fitted amplitudes, which gives an idea of how accurate our fit is. 

For the honeycomb lattice high-temperature series, we conduct two separate fits, one at 
the ferromagnetic point (with additional suppression of the antiferromagnetic singularity) 
and one at the antiferromagnetic point (with additional suppression of the ferromagnetic 
singularity). In fact we also did this for the square lattice, to check for an antiferromagnetic 
scaling term. We found no such scaling term, which is consistent with the results in 
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Figure 5: The conjectured natural boundary in the complex (cube root) nome g-plane for 
the Ising model on the triangular lattice. The real axis cusps are the points u = ±1/3; 
the other two are u = —1 + iO^. The right side arcs are defined by |A:tr| = 1 and u = 
(— 1 + 2e"'^)/3, — TT < 4> < TT. The left arcs correspond to straight line segments lying on 
either side of the cuts, — oo < fctr < —1 in fctr and — 1 < u < — 1/3 in u. Crosses mark the 
singularities found in the series analysis described in the text. For clarity, the singularities 
for T < Tc are shown only in the upper half plane, those for T > Tc in the lower half. 
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Once the initial fitting has been done, we can improve the accuracy of our fits by 
iteratively subtracting the new fit (or fits), re-suppressing singularities (replicating this in 
our fitting functions) and fitting again to the remainder, and so on. 
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Appendices 

A Ferromagnetic scaling function 
A.l Squaire lattice 



Fl^ = (t + Vl + r2)^(l + TV2-rVl2 

- 6.3213068404959366230670987124576163379333404464\ 

29429335850509012099708742399 • 
+ 6.2519974704602432856837331806319562265626657486\ 

9581059930911004970341 • 

- 5.6896599756179940495694760341390552949459234168\ 

0072164185003897 • r^" 
+ 5.14221827114214604273511179366558788399868131986546472359 • r^^ 

- 4.67471611538219753943422533513538091798878146367647 • r^^ 
+ 4.28351401741664147913747092020949150840022385 ■ r^^ 

- 3.93463085065515612248985707350481524149 • r^^ 
+ 3.613033718221972872129117995447426 • 

- 3.3030941616500642890625665822 • r^^ 

+ 2.99419136711436481655789 • r^^ - 2.674815242128336541 • r^^ 
+ 2.3339198769874 • r^^ - 1.95837351 • + 1.537 • r^^), 



= (r + v^rT72)^(l + rV2-rVl2 

- 0.12352922857520866639356466570562347322323268198504142433416176 • r' 
+ 0.13661094980909643478343857458083310826834711524701276519 • 

- 0.13043897213329076084013583556244683622929916938362 • 
+ 0.121512875791442694842447521021056149318718395 ■ r^^ 

- 0.1129603634344171840043033744010408654148 • r^^ 
+ 0.10536961142693738687373469324338873 • r^^ 

- 0.0982140320131209895954107399728 • r^^ 

+ 0.091314688764698386593329786 • r^^ - 0.08439419183682814997218 • r^^ 
+ 0.0772604004964458205 • r^^ - 0.069668638313388 • r^^ 
+ 0.061368727265 ■ r^^ - 0.05204288 ■ + 0.0414 • t^^). 
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Triangulcir lattice 



F^J = (r + + r2)5(l + 1/2 . - 21/256 • 

- 6.7764559898170749532861771919188746477857219070(3) • 
+ 6.84262914118601551543582352085826620764414(10) • 

- 6.250933162702506214104998011755062095(9) • 
+ 5.63987692190321788346983658716286(30) • r^^ 

- 5.106253322544511659092052061(5) -T^^ 
+ 4.65493974449161799368079(6) • r^^ 

- 4.2701171199002454178(4) • r^^ 

+ 3.9327480363388237(23) • - 3.625158242566(11) • r^^ 
+ 3.33306138(7) • r^"^ - 3.04765(11) • r^^), 



Fj"^ = (r + Vl + t2) i (1 + 1/2 • - 21/256 • 

- 0.1359799770448664282788192846845965785(4) • 
+ 0.152349558318015426490910429319733(17) • 

- 0.14450411683821267150571729255(18) • 
+ 0.1331875171226390774852445(8) • r^^ 

- 0.1223854265244510620558(16) • r^'' 
+ 0.1128620837499335229(18) • r^^ 

- 0.1045232876841806(12) • r^^ 

+ 0.0970484952533(5) • r^^ - 0.09008180554(18) • r^^ 
+ 0.08331757(8) • r^^ - 0.07654(4) • r^^). 
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Honeycomb lattice 



F^^ = (r + Vl + + 1/2 • - 21/256 • 

- 2.2311493924390249844287257306396248825952406357(14) • 
+ 2.29732254380796554657837205957901644245366(35) • 

- 2.169834272110400332644049880573751163(24) • 
+ 2.0232401262820557301956317745273(7) • r^^ 

- 1.887361520807880372000774531(10) • r^^ 
+ 1.76703614250551894208334(10) • r^^ 

- 1.6608531413942897306(6) • r^^ 

+ 1.5669217708308492(27) • - 1.482989258248(11) • r^^ 
+ 1.40629074(6) • r^^ - 1.33479(6) • r^^), 



= (r + Vl + r2)§(l + 1/2 • - 21/256 • 

- 0.01765738818162214275960642822819883(11) ■ 
+ 0.0340269694547711409716975728625(27) • 

- 0.038525626277127202618271411(17) • 
+ 0.03947785106181932194262(4) • r^^ 

- 0.03910559858848358918(5) • r^^ 
+ 0.03820301383607213(3) • r^^ 

- 0.037081202839025(16) • r^* 

+ 0.035888270323(6) • r^" - 0.0346847083(13) • r^^ 
+ 0.03347384(20) • r^^ - 0.032213(27) • r^^). 
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B Antiferromagnetic scaling function 
B.l Honeycomb lattice 

^hcjaf ^ -(r+ Vl + r2)5 

X (4.545306597378049968857451461279249765190481271258(18) • 

- 4.545306597378049968857451461279249765190481271258(18) • 
+ 4.0810988905921058814609481311813109325(5) • 

- 3.61663679562116215327420481263559(4) • r^^ 
+ 3.2188918017366312870912775304(12) • r^^ 

- 2.887903601986099051597450(19) • r^^ 

+ 2.60926397850595568720(21) • r^* - 2.3658262655079745(17) • r^" 
+ 2.142168984318(13) • r^^ - 1.92677064(12) • r^^ + 1.7124(8) • r^^), 



^hcjaf ^ -(r + v/1 + t2)5 

X (0.1183225888632442855192128564563977189(6) • 

- 0.1183225888632442855192128564563977189(6) • 
+ 0.10597849056108546888744587531(10) • 

- 0.0937096660608197555426090(10) • r^^ 
+ 0.083279827935967472857(4) • t^'^ 

- 0.074659069913861378(7) • r^^ + 0.067442084845148(6) • t^^ 

- 0.061160224928(3) • r^^ + 0.0553970966(11) • r^^ 

- 0.0498435(28) • r^^ + 0.04424(7) • r^^). 

Comparison of these results with what is obtained from the ferromagnetic expressions 
in Appendices A. 2 and A. 3 using yields a partial check of the consistency of our 



numerical fitting described in Section |4j 



C Short-distance terms 

Here we give the short-distance "regular" background terms of the form 

oo Lv^J 

rounded to 15 places. Our complete results are available in [33]. 
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C.l Ferromagnetic square lattice 



Ssq = (r + yrT7^)V2 

X [- .104133245093831 - .074368869753207 r - .008144713909120 
+ .004504107712232 + .239618794254722 r'' - .002539950595339 

- .235288909669962 + .001915707531701 + .214340096611538 

- .000883215706003 - . 194220628407196 + .000007233509777 
+ . 177102037555467 T^^ + .000688811096268 r^^ - .162792536489746 

- .001236572355315 r^^ + .150013412064378 r^^ + .001671694059110 r^^ 

- .138208109106217 r^^ - .002022002972782x1'' + .126799277310505x^0 
+ . 002308285588780 r^i - .115396441906289x^2 - .002545765264414x^3 
+ .103574086263807x^4 + .002745532102527x^5 - .090922989554413x2^ 

- .002916073299270x2^ + .076954225263348x2^ + .003063568441388x2^ 
+(ln|x|) 

X (+ .032352268477309 X - .005775529379688 x^ + .059074961290345x4 

+ . 003058491575856 x'5 - .059166272208841 x^ - .002067088393167 x'' 
+ . 054246930704214 x« + .001060102531550 x'' - .049300253157083x1° 

- .000268300641612 x" + .045027052571957x^2 - .000343326832572x^3 

- .041428586463053 x" + .000819393297118 x^^ + .038202673904453x1^ 
- .001196464684146x1^ - .035217475800642x1^ + .001500680711946x1^ 
+ .032331741680806x20 - .001750700134389x21 - .029449221445927x22 
+ .001959866123653x23 + .026464090269923x24 - .002137779981361 x25 

- .023274239921560x2*5 + .002291702790868x2^ + .019757464449312x2^ 
-.002426942629382x2^) 

+(ln|x|)2 

X (+.009391569871146x4 - .008695925462879 x*^ + .007669481493105 x^ 
+ . 000154284382979 x^ - .006805407688144x1° - .000310520937481x11 
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+ . 006113866432195 r^^ + .000444606198236 - .005557100215116 

- .000554418149346 r^^ + .005078042485427 r^^ + .000643607994970 r^^ 

- .004649202184071 r^* - .000716232782651 r^^ + .004246382079429 r^" 

+ . 000775832889819 - .003853404958387 r^^ - .000825213786325x23 
+ .003454329481031x24 + .000866512510954x25 - .003034537504706x26 

- .000901440561715x2^ + .002577451310655 x2« + .000931250046525 x23) 
+(ln|x|f 

X (- .000015771569138 x^ + .000034428206621 x" - .000052442717749x^3 
+ .000068823835730x^5 - .000002084325090x^6 - .000083482363640x1^ 
+ .000006458964601x1* + .000096589603855x1^ _ .000013639281329x2" 

- .000108385585447x21 + .000023853448397x22 + .000119105615864x23 
- .000037600547029x24 - .000128947973257 x25 + .000055460969100 x26 
+ .000138099034068x2^ - .000078321412692x2^ - .000146701272364 x29) 

+(ln|x|)* 

X (- .000000145427323x16 + .000000452982068x1* - .000000959267146x2° 
+ .000001683186013x22 - .000002660926741 x2'i - .000000003368087 x25 
+ .000003934622630x26 + .000000009693809x2^ _ .000005565949306x2* 
-. 000000023894457 x2'J) 
+(ln|xi)5 

X (+ .000000000141953x25 - .000000000441519x2^ + .000000001224727 x29)]. 



C.2 Antiferromagnetic square lattice 

B^l = (r+^/l^)i/2 

X [+. 158866522960947 +. 149566836938536 X + .010712225879833x2 
+ .012753018839962x3 - .011741188869656 x^ - .014066040875666x5 
+ .013106454615626 x6 + .012239696625538 x^ - .011840194045411 x* 

- .010585409302312 x^ + .010151560037724x1" + .009080004112331 x" 

- .008542012228790x12 - .007717026940132x13 + .007123677682511x1' 
+ .006508646391366x15 - .005912245104109 xi^ - .005467010793273x1' 
+ .004900133679335x1* + .004597213578999x1^ - .004074131287647x2' 

- .003893839411793x21 + .003417128190380x22 + .003339548120697x2' 

- .002905973440848x24 - .002906261172134x25 + .002510579952576x2' 
+ .002559795034096x2^ - .002197017191525x2* - .002268131101616x2' 



43 



+(ln|r|) 

X (- .155317190158011 r + .032067148145870 - .007716887572462 r'' 

- .015675211573817 - .000285542451537 r*^ + .009607254502732 

+ . 004835406420625 - .006064990344481 - .007340015041447 r^" 
+ . 003910356521404 + .008708427445682x^2 _ .002697783010885 r^-' 

- .009405623038077x1'' + .002161267525775 r^^ + .009683424894714 r"' 

- .002106502189570x1^ - .009696425760611 r^^ + .002374179655585 r^^ 
+ . 009556527066075 - .002827702235523x^1 - .009355445823856x^2 
+ .003353500098021x^3 + .009167728425385x^4 - .003868501274293 x=^^ 

- .009041887308405x^6 + .004329857870148x2^ + .008988699114041x2^ 
-.004740651061453x29) 

+(ln|x|)2 

X (+ .011533714378823 X* - .011311734920692x6 + .010045768711199 x^ 

- .000475698571097 x^ - .008783972022287x1° + .001157180172964x11 

+ .007680651109513x12 - .001865091261620x13 - .006744701894515 x" 
+ .002491836298308x1'^ + .005964068368078 xi^ - .002972695839442x1^ 

- .005327647984236x1^ + .003273731192324 xi^ + .004822634345908x2" 

- .003388237983845 x21 - .004425467442049 x22 + .003335544654325 x23 
+ .004094631068058x24 - .003157926493823 x25 - .003772517361799 x26 
+ .002912409430061x2^ + .003402363586667x2^ - .002655520057041 x29) 

+(ln|x|)3 

X (+.000057899719476x9 - .000169915088240 r" + .000326648846875 xi^ 

.000517190858645x1^ - .000001422188017 xi^ + .000729027463661 x" 

+ .000009170599968x1^ - .000948102082432 xi^ - .000032108584334x2° 
+ .001159641637018x21 + .000080604994072x22 - .001349862803424 x23 

- .000160502424863x24 + .001508463556122 x25 + .000264448046627 x26 
- .001631439353298x2^ - .000364294406680 x28 + .001722930696448 x29) 

+(ln|x|)4 

X (- .000000160856746x16 + .000000456983407x1^ - .000000040918655x2° 

- .000003417322708x22 + .000013997416073 x24 - .000000009021984 x25 

- .000036813182410x26 + .000000125169178x2^ + .000075735555538 x2« 
-.000000779435552x29) 

+(ln|x|)5 

X (- .000000001286222x25 + .000000004385050x2^ + .000000027799861x29)]. 
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C.3 Ferromagnetic triangular lattice 



St, = (r+ Vi + t2)1/2 

X [- .049561116521763 - .029358763163227 r - .003802085786368 
+ .006390376143904 + .194331491416170 r"* - .004659659320547 

- . 195488838278358 + .003651173504528 + .178621656686715 

- .002895748949957 - . 161336242614720 r^" + .002311455809247 r" 
+ .146284971711630x^2 _ .001842503338574x^3 - .133577942096403 

+ .001456977882233x^5 + .122753362974429 r^^ - .001134805649081 r^^ 

- .113263489337451 r^*^ + .000862240813452 r^^ + .104601807098273x^0 

- .000629270080307x^1 - .096359884476827x^2 + .000428320304385x^3 
+ (ln|x|) 

X (- .005374288589598 X + .001021325616916x3 + .049253501657254 x"* 

- .000006005387528x5 - .050675128993180 x^ - .000277768605459 x^ 
+ .046680337431830 x^ + .000300252836069 x^ - .042334826787302x1° 

- .000227268452048 x" + .038476859060257x^2 + .000114673447726 x" 

- .035187926212720 x" + .000011949942590x^5 + .032370727904288x1^ 

- .000140047340708 x" - .029894415736149x1^ + .000263535279625x1^ 
+ .027634040948919x2" - .000379722086737x21 - .025487516010638x22 
+ .000487692500440x23) 

+ (ln|x|)2 

X (+ .008301571737990x4 - .007863822472801x6 + .006940825976817 x^ 
+ .000004920887586 x^ - .006124967722414x1° - .000008028657674 x" 
+ .005459215424842x12 + .000010897860945x13 - .004918734820641 x" 
- .000014290481783x15 + .004471830270400 xi^ + .000018331038938x1^ 

- .004090720207698x1** - .000022918919968 xi^ + .003752754837032x2° 
+ .000027921194209 x21 - .003440642618124 x22 - .000033228452894 x23) 

+ (ln|x|)3 

X (- .000008243826432x9 + .000019339918959 x" - .000030791103742x13 
+ .000041558105858x15 - .000000289107374 xi^ - .000051264920936x1^ 
+ .000001306093325x1** + .000059839322918 xi^ - .000003538519413x2° 

- .000067344361243 x21 + .000007394079439 x22 + .000073895527830 x23) 
+ (ln|x|)4 

X (- .000000023041822x16 + .000000103665399x1* - .000000280140679x2° 
+ .000000584762202x22)]. 

The leading term in iJtr above confirms the estimate in (26) of [2D] after adding a factor 2 needed 
because of a difference in conventions. 
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C.4 Ferromagnetic honeycomb lattice 



She = (T+v/l + r2)V2 

X [- .221526277068482 - .170518806873542 r - .019236029093417 

- .000240258087320 + .140112831065240 + .002715126912573 

- .139853710977321 - .002422123613147 + .130966255841349 

+ .002478441994620 - .121767984156181 - .002593394970413 r^^ 
+ .113507041635338x^2 _^ .002697199234761 r^^ - .106296857456039 

- .002778848577189 T^-^ + .099986481192228x1*' + .002840749390802 

- .094433287150951 r^*^ - .002887196461298 r^^ + .089489266542983x2° 
+ .002921915790808 - .084984336381365x^2 - .002947776518703x^3 

+(ln|T|) 

X (+ .110304596706594 T - .017367191250168x3 + .032554394731493 
+ . 007749610093406 - .033370773266168 r'^ - .004545306065368 

+ .031517394252347 + .002697183732340x9 - .029438657077664x1" 

- .001512880972029 x" + .027527821013048x12 + .000704802950233x13 
- .025841562449391 x" - .000126907282155x1^ + .024355378143637x1^ 

- .000302184134792x1^ - .023041033316250x1® + .000630617459438 xi^ 
+ .021866143676416x2° - .000888340903130x21 - .020791271125817x22 
+ .001094822052575x23) 

+(ln|x|)2 

X (+ .004328421950579 x'' - .004173174221495 x*' + .003864856815018 x® 
+ .000098594831882x9 - .003581293629898x1° - .000203109598421 x" 
+ .003336710555165x12 + .000293321457772x13 - .003126961910309 x" 

- .000367344803112x1^ + .002944831872355 xi'^ + .000427433155150x1^ 

- .002785081058367x1® - .000476286169756 xi^ + .002642924163304x2° 
+ .000516221811618x21 - .002513107496440x22 - .000549081268363x23) 

+(ln|x|)3 

X (- .000008459084030 x° + .000018328291997x11 - .000027731894823x13 
+ .000036311733546x1^ - .000000100196269x1^ - .000044084313016x1^ 
+ .000000291941736x1® + .000051136242275 xi^ - .000000616104491x2° 
- .000057556878974x21 + .000001226160802x22 + .000063427686583x23) 

+(ln|x|)4 

X (- .000000008316207x1° + .000000024195538x1® - .000000049128873x2° 
+ .000000092415863 x22)]. 
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C.5 Antiferromagnetic honeycomb lattice 



B^^ = {t+ \/l + T^)^/^ [.122404044024957+ .1118012805470871 
+ (ln |t|) (- .121053173885789 r + ...) + ...], 



as given more fully by (13) and Appendices C.3 and C.4 
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